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INTRODUCTION 


A three  dimensional  (3D)  finite  difference  code  is  be- 
ing developed  which  will  perform  a deterministic  simulation 
of  a stick  slip  earthquake.  The  primary  goal  of  this  research 
effort  is  to  obtain  an  earthquake's  equivalent  elastic  source 
in  a form  which  is  suitable  for  propagation  to  teleseismic  dis- 
tances. Normalization  of  the  earthquake  model  to  the  near 
field  data  obtained  at  Bear  Valley  is  expected  to  furnish  the 
required  validation  of  the  model. 

RkSULTS 

In  order  to  obtain  a 3D  code  that  is  both  computationally 
efficient  and  still  flexible  enough  to  permit  nonlinear  material 
response  in  the  fault  -one,  a technique  was  developed  which 
links  a special  purpose  small  deformation  linear  elastic  code 
to  a nonlinear  material  response  code. 

This  link  was  tested  in  two  dimensions  by  merging  CRANl/2^ 
which  contains  the  two  dimensional  fault  model  reported  by 
Cherry, ^ and  a special  purpose  linear  code  which  assumes 
a Hooke's  Law  Material  behavior.  The  difference  equation  used 
in  CRAM  to  simulate  the  fault  surface  and  the  assumed  nonlinear 
material  behavior  in  the  fault  zone  are  given  in  Appendices  A, 

B and  C.  The  difference  equations  for  the  2D  linear  code  are 
given  in  Appendix  D. 

Calculation  5A  was  used  as  the  test  problem.  This  cal- 
culation1 ^ featured  a 5 km  fault,  a rupture  velocity  of  2.15 
km/sec  and  a dynamic  stress  drop  of  0.5  kbar.  A schematic 
of  the  calculation  is  shown  in  Fig.  1.  This  calculation  was 
rerun  using  the  linear  code  linked  to  CRAM.  The  CRAM  grid 
was  a rectangle  with  dimensions  13  km  x 7 km,  sufficient  to 
>.over  the  fault  and  the  nonlinear  region  in  the  fault  zone. 
Calculations  outside  this  rectangular  region  were  performed 
by  the  linear  code. 

Comparison  between  the  original  calculation  5A  and  the 
new  calculation  in  which  the  linear  code  was  linked  to  CRAM 


1 


is  given  in  Figs.  2 through  11.  Agreement  between  the  two 
calculations  is  excellent. 

The  motivation  for  using  the  linear  code  for  calcula- 
tions in  the  small  displacement  elastic  regime  is  to  reduce 
both  computation  time  and  computer  costs.  The  new  calcula- 
tion, using  the  link,  reduced  the  cost  of  the  simulation  by 
a factor  of  4.3. 

A three  dimensional  (3D)  linear  elastic  code  was 
developed  (Appendix  E)  and  the  framework  for  linking  to  a 
3D  nonlinear  code,  containing  the  stick-slip  fault  model, 
was  included  in  the  linear  code. 

In  order  to  debug  the  3D  code  comparable  problems  were 
run  on  the  2D  and  3D  codes.  Figure  12  shows  the  regions  over 
which  surface  tractions  were  applied  in  the  two  codes.  The 
21)  code  was  run  in  axisymmetric  geometry  and  the  surface  trac- 
tion region  in  the  3D  code  was  made  as  nearly  circular  as 
possible . 

Figure  13  compares  the  radial  component  of  particle 
velocity  between  the  two  codes  at  50  usee  while  Fig.  14  com- 
pares the  vertical  component  of  particle  velocity.  Differences 
near  the  source  are  probably  caused  by  the  approximations  to 
the  circular  area  of  loading  required  in  the  3D  simulation. 

In  its  present  form  the  3D  linear  code  is  capable  of 
determining  the  modification  of  the  free  field  earthquake 
ground  motion  produced  by  a non-homogeneous  geologic  environ- 
ment. A given  equivalent  sourep  may  be  used  to  drive  the  3D 
code  in  order  to  obtain  theoretical  seismograms  at  seismo- 
meter locations  in  the  near  field.  Linking  of  the  3D  linear 
code  to  a nonlinear  code  containing  a stick  slip  rupture 
model  is  proceeding  under  separate  funding. 
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Figure  2.  Particle  velocity  at  Station  1,  calculation  5A 
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Figure  5.  Particle  velocity  at  Station  1,  calculation  5A  re- 
peated with  linear  code  linked  to  CRAM. 
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Figure  4.  Particle  velocity  at  Station  2,  Calculation  5A. 
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Figure  10.  Particle  velocity  at  Station  5,  Calculation  5A. 
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Figure  12.  Surface  areas  over  which  the  pressure  load  was 
applied  for  2D  and  3D  comparison  calculations. 
The  cross  hatched  area  corresponds  to  the  3D 
calculation. 
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APPENDIX  A 

The  difference  equations  used  in  CRAM1  1 to  move  an 
interior  point  are  written  such  that  the  boundary  conditions 
for  an  exterior  point  are  obscured.  Since  we  would  like  to 
isolate  the  normal  and  tangential  stresses  at  an  interior 
interface,  this  requires  that  the  interface  be  treated  as  an 
exterior  line  over  which  the  correct  boundary  stresses  are 
applied.  In  this  section  a differencing  scheme  is  obtained 
that  isolates  the  boundary  stresses  and  that  is  consistent 
with  the  CRAM  interior  difference  equations. 

The  conservation  of  linear  momentum  (equation  of  motion) 
in  two-dimensional  Cartesian  geometry  is 


If  a Lagrangian  coordinate  system  (k,j)  is  established  in 
the  material,  then 


9Z  _ 3Z  9x  9Z  3y 

9l<  9x  Tl<  Ty  Hk 


91  9E  9x  9E  9y 

9j  7x  TJ  Ty  9j 


where  E is  a typical  stress  component  (xx,  yy,  xy)  in  the 
equation  of  motion.  Solving  for  9E/9x  and  9E/9y  gives 


9 E _ 1 ["  9 E 9y  9 Z 3y ' 

cTx  J 9j  3l  TF  TjT 


9 Z _ _ 1 9 E 9x  _ c E 9x 
Ty  T 9j  9T  9j 


(3) 


1 s 


(4) 
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where 


T _ 3x  3y  3x  3v 

J " 3J  TF  " Tic  Tj 

If  the  k,j  coordinates  assume  discrete  values 


then 


u here 


1,  2,  ...  k- 1 , k,  K+l,  ...  k 

max 


I>  *-»  •••  j*l»  j,  j+l> 


\ = | j Rj.  > ina  e = 2Aa  e 


[5  = A + y -A  t - 3xo  .9v-> 

Rj  Tf  Cx  3j  ej  V 31  ex  + It  ej 


Also 


i$.  X 5,  = (1*  . 3x  3^\  - + 

J k V 3 j 37  TfT  3j/  ex  ey  J 


e,  x e 
k y 


Comparing  Eqs.  (6)  and  (7)  shows  that  a good  approximation 
to  J will  be  the  zone  area  (Aa  + a^)  if 


e = e x e 
x y 

Equation  f 8 ) is  satisfied  if  the  x,  y and  k,j  coordinates 
have  the  same  relative  orientation  as  that  shown  in  Fig.  1, 
i.e.,  if  the  unit  vector  obtained  9 8 7 

from  the  j » I operation  is  equal  / T^s^b 

to  the  unit  vector  from  e^  * e^.  b 


y 


* i 


Fig.  1--The  x,y  and  k,j 
coordinate  system. 
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i-igure  i also  shows  tne  nuniDering  system  to 
oth  nodal  point  and  interior  variables.  Table 
quivalence  between  the  numbering  system  and  th 


Number 


k,j  Interior 


k-1/2,  j - 1/2 


k-1/2,  j +1/2 
k+1/2,  j +1/2 
k+1/2,  j+1/2 


k , j Exterior 


k,  -1 
k-1,  j-1 
k-1,  j 
k-1,  j + 1 
k,  j+1 
k+1,  j+1 
k+1,  j 
k+1,  j-1 


Typical  interior  variables  are  density,  stress, 
internal  energy  and  area.  Exterior  variables  ar 
tion,  velocity  and  position  vector. 

Equations  (3)  and  (4)  and  Fig.  1 si 
natural  differencing  scheme;  i.  e., 


y 

7 8 8 1 


P J + (1 
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I i I 

p 


£ X 
7 8 8 1 

P J + P J 


W, 


£ x 

6 I I 4 


I X w . 


7 7 


8 8 


where  £ = £ 

2 8 7 


8 8 1 


pT 
6 8 

+ 

p 

i i 

~ 7 

£ 

X 

8 

i 

1 2 

P J 

+ 

P 

8 8 

i i 

T 

= >\ 

_ 

y 

» 1 

8 

i 

p~3  + p J 

7 7 6 6 

2 


(1 


V 


(10) 


- • y C1J1U  w CUIU  W • 

weight  the  individual  acceleration  components  based  on  thei 
location  with  respect  to  point  1. 

A fairly  simple  weighting  scheme  used  in  the  TENSOR 


code 


[3] 


is 


w. 


a • a 

14  8 4 


a • a 

8 4 8 4 


Uj  " a 


a • a 

12  6 2 


6 2 


6 2 


In  Eqs . (9)  and  (10] , if 


Kk  = wj  = 1/2 


(id 


Pcjc  + Pdjd  P-J.  + P-J-  + P-J-  + p J 


T 


d d 


i i 


2 2 


3 3 


4 4 


(12) 

[4] 


then  the  CRAM  interior  differencing  is  obtained.  Both  HEMP 
and  CRAM  use  the  same  interior  differencing  scheme.  Equations 
(11)  and  (12)  reduce  the  TENSOR  difference  equations  to  those 
of  CRAM  and  HEMP. 


If  a k line  is  to  be  decoupled  from  the  grid,  as 
in  Fig.  2,  then  Eqs.  (9)  and  (10)  permit  this  if  = 1 
for  k and  w^  = 0 for  k . For  a point  on  k+,  w^  = 1 
and  vjz  = 1/2.  Equations  (9),  (10)  and  (12)  may  be  used  to 
write  the  spatial  derivatives  in  Eqs.  (1)  and  (2)  giving 
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Figure  ...  A decoupled  interior  grid  line.  The  boundary 
stresses  kk*  and  Tcj”*  do  not  change  when  the  interface 
is  viewed  from  below  (w-j.  = 1)  or  above  (wk  = 0). 
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Figure  3.  The  orthogonal  j,  k unit  vectors. 
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1 3 xx 

p 17 


1 3xv 
P 3y 


xx  y 

7 8 8 1 


xy  x 


(XK,  ' 


xx*  y + | xx  - xx 
£ ' 6 1 \ 8 


* )y 

1 / 12 


7 8 8 1 


I 

p 3y 


yy  * 


7 8 8 1 


I xy  - xy*  )x  + I xy  - xy*\  x 

) 7 6 1 6 1 \ __8 1 / 1 

(77,  - yi>)xti  * (77,  - Ty*)*, 


l 3x2  . !^Za. 

p 3x  + 


( xy  - xy*)y  + ( xy  - xy*  )y 

1 7 6 ' 6 1 \ 8 1 / 1 ? 


where 


P J P 

+ _ 7 7 


J 

8 8 


Along  a typical  interface  line  joining  two  adjacent 
nodal  points  (Fig.  3)  orthogonal  unit  vectors  £ and  j 


k = -ev  sine  + e cos0 
x y 


J = e cose  + e sine 
a y 


The  stress  components  in  this  coordinate  system  are 


Fk"  = yy  cos20  + xx  sin20  - 2xy  sine  cos6 


FJ  = xy(cos26  - sin2e)  + (yy  - xx)  sin6  cose 


jj  = xx  cosz6  + yy  sin26  + 2xy  sine  cosO 


The  acceleration  of  a point  on  k may  be  written 


'dx\+  1 

\Ji!  = i 

<t> 


7 8 6 


6 6 1 


xx  y + xx  y - xy  x - xy  x + FIc*y+ 

828  786  828  * * 

]<F*y*.-  FJJx^-  kj*x"2] 


(13) 

(14) 

(15) 

(16) 

(17) 
are 

(18) 

(19) 

(20) 
(21) 
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(22) 
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dv  + 

It 


-yy  x 

7 8 6 


yy  x + xy  y + xy  y - Ec*x 

828  786  828  661 


- kl*x+  - Tc7*y+  - TcT*y+ 

1 12  661  112 


(23) 


Similarly,  the  acceleration  of  a point  on  k"  may  be  written 

as 


dx  \ 1 

It)  = “t 


XX  y + XX  y - xy  x - xy  x - Ec*y‘ 

6 6 4 1 ‘t  2 664  142  661 

- kT*y  + FjT  x’  + l<7*x~  1 

1 1 2 6 6 1 i 1 2 J 


(24) 


dy  = i 

' A 


It 


4> 


-yy  x - yy  x + xy  y + xy  y + Ec*x‘ 

664  142  664  1 k 2 661 

+ kk’*x  + lcj*y  + Ff*y"  1 

1 12  6 6 1 1 1 2 J 


(25) 


These  equations  have  been  derived  assuming  that  Tele*, 
kj*  kk*  and  kj*  are  specified  along  the  6-2  interface, 

oil  7 

finding  the  corresponding  stresses  in  the  x,y  coordinate 


system,  i.e. , 

xx*  = yy*  cos2e  + £7* 

sin2  0 

- 2l7*sin0  cos0 

(26) 

yy*  = 37*  sin20  + ktc* 

c o s 2 0 

+ 2FJVisin0  cos0 

(27) 

xy*  = (37*  * ^*)sin0 

cos  0 + 

3T*(cos20  - sin20) 

(28) 

and  then  substituting  these  expressions  for  xx*,  yy*,  xy*, 
xx*,  yy*,  xy*  into  Eqs.  (13)  through  (16) 

Equations  (22)  through  (25)  will  be  used  to  move  the 
points  on  the  decoupled  grid  line.  They  are  consistent  with 
the  interior  difference  equations. 
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APPENDIX  B 


In  order  to  use  Equations  22  through  25  in  Appendix  A 
the  boundary  stresses  must  be  specified.  In  order  to  solve 
for  these  stresses  we  apply  contact  discontinuity  boundary 
conditions,  which  require  that  the  normal  component  of  stress 
and  normal  component  of  velocity  be  continuous  at  the  bound- 
ary . 


The  unit  vectors  normal  and  tangent  to  the  6-2  inter 
face  may  be  written 


le  = - s in0  c + cose  e 
x y 


j = cose  e„  + sine  e 
x y 


e^  = -sine  t + cose  J 


e^  = cose  le  + sine  J 


CD 


The  acceleration  components  on  the  plus  side  of  the  line  may 
be  written 

ak  = ^t[2l  + R ^ + R Be  - (R  +R  ) Bc*l 

k C+)  L k 6 1 7 1 2 6 6 1 1 2 J 

ai  = ^rfg]  + R + r rj  - (r  +r  ) ft*| 

J (j,  LJ  61  7 12  e 61  12  J 

The  corresponding  acceleration  components  on  the  minus  side 
are 

a"  = — [g.'  *•  (R  +R  ) Be*  - R Be  - R Be] 
k 4)'  L6k  6 112  6 16  1 2 1 J 

a:  = ~[g-  + (R  +R  ) EJ*-  R Bf  - R lg  ] 

1 (+)Lj  6112  616  12  1 J 


(2) 


(3) 


(4) 


(5) 
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Equations  2 through  5 have  been  derived  from  Equations 
13  through  16  in  Appendix  A by  assuming  that  the  boundary 
stress  is  uniform  over  the  entire  6-2  portion  of  the  grid 
line.  These  equations,  therefore,  are  not  completely  con- 
sistent with  the  interior  difference  equations.  The  assump- 
tion concerning  uniform  boundary  stress  is  necessary  since 
the  contact  discontinuity  boundary  conditions  will  result  in 
an  equation  that  relates  the  normal  components  of  acceleration 
on  the  plus  and  minus  side  of  the  boundary.  This  equation 
should  contain  only  one  unknown,  i.e.,  the  normal  component 
of  the  boundary  stress. 

In  order  to  find  the  relation  between  the  above  acceler- 
ation components  at  a slipping  interface,  we  follow  the  tech- 
nique used  by  Cherry,  e t a 1 , ^ ^ and  attach  a coordinate  system 
to  the  point  on  the  minus  side,  with  k and  J being  the  unit 
vectors  normal  and  tangent  to  the  slip  line  at  the  point  to 
be  moved. 
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If  v and  v are  the  velocities  on  the  plus  and  minus 
side  of  the  slip  line  then,  from  Equation  1,  the  normal  com- 
ponents of  velocity  are  given  by 

■+■+  :►  • + . • + 

v • k = -x  sin9  + y cosG 

(9) 

T*  • - • - 

v • k = -x  sinQ  + y cosG 
If 


V ' = 9x  ' £Z. 

Z dZ  YZ  = 3 Z 


then 


y'i 

Sin0  = rr—  COS  0 

z 


(10) 


Substituting  Equation  10  into  9,  and  equating  normal  velocity 
components  gives 

= V+-k  (id 


or 


= -*\*y*xi  az) 

Since  the  k,l  coordinate  system  is  attached  to  the  minus  side 
of  the  slip  line  then 


d 

dt 


C-x'y^y'x") 


3 . • + - • + - 

■ n t_x  V* 


k , fi.  constant 


(13) 
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Ihe  left  side  of  Equation  13  may  be  written 

•>'?,  dtCx  ) + x£  dt y * x x£ 

= R^ak  ♦ (v  * ^3  Cv£ * J ) - (v'*I)(v'*£) 


(14) 


Since 


- !r(v+  -j)  ♦ V ^ 


I_r£  y.t) 

R 3k  V M 


'k 


and 


(v  + - v)  • k = 0 


Ihcn  the  right  side  of  Equation  13  may  be  written 


• + - • + - . 


atc'x  y*+>’  x£) 

= R£ak  + (v+*l<)  ( y • J)  - (v'.I)(v'*£)  - (v+-v‘)  *£v+*i< 


(15) 


In  Equations  14  and  15  a^  and  a^  are  given  by 


a*  = - }Jl  + 1± 

k 3t  R^  dt  Rz 


a;  = - j-tiO  — + h. 

k dt  R^  3t  R^ 
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and  arc  the  sane  acceleration  conponents  given  by  Equations  2 
and  4.  We  have  also  used  the  relation 


-X  «, 


(v--kKv£*j) 


(v+  -t) 


C16) 


in  order  to  derive  Equations  14  and  15. 

Substituting  Equations  14  and  15  into  Equation  15 

gives 

+ 

a,  -a,  = A 
k k c 

who  re 


(16) 


A 

C 


(V  -V  ) -l (v5+vc ) • k 


(17) 


Solving  Equation  16  for  kF*  gives 


FF*  = 


0 g^.-b^g,  + $’(R  Fk +R  FF  ) 

k K 6 1 7 12  6 

(R  +R  )($++0'j 
6 1 12 


q+  (R  Fk +r  FF  ) - d+0  A 
6 1 6 12  1 ¥ c 

(18) 


Equation  18  relates  the  normal  component  of  stress  (Fie*) 
at  the  slipping  interface  to  interior  zone  variables.  In  order 
tc  move  the  boundary  points  then  FJ*  in  Equations  22  through  24 
in  Appendix  A must  also  be  specified. 


For  a tied  point  A in  Equation  18  is  set  equal  to  zero, 

s l nee 


(vT  -\T)  • £ = 0 

Also  the  tangential  stress  component,  for  a tied  point,  is 
obtained  from  Equations  3 and  5.  Since 
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a a 

J J 

then  for  a tied  point 


+ 0 (R  FJ  FJ  ) + o+(R  FJ +R  FJ  ) 

J J 617  12  B 6 1 _ 6 12  J l 

(R  +R  ) (4>++o ' ) 

6 1 12 


FJ*  = 


K- £4  7 o 


APPl'.NDl  X C 


Plastic  flow  is  due  to  the  inability  of  real  materials 
to  support  unlimited  values  of  shear  stress.  In  the  code  the 
deviatoric  stress  components  are  modified  such  that  the  re- 
sulting stress  state  is  consistent  with  a Mises  yield 
criterion. 

If  the  second  deviatoric  invariant  (J)  is  greater  than 
a specified  value  (1/3  Y2) , then 


s..  - §..  JL  (j  > vi) 
1J  1J  yjj  \ > ) 


(1) 


where  S - j is  the  adjusted  stress  deviator 

A 

S—  is  the  stress  deviator  calculated  by  assuming 
that  the  total  strain  rate  is  elastic,  and 


j ■ ; csjj  sjj) 


(2) 


For  a triaxial  test,  Y corresponds  to  the  maximum  allowable 
stress  difference  at  failure. 


Rupture  initiation  is  being  modeled  by  accumulating 
the  difference  between  /IT  and  Y during  yielding.  When 
this  accumulation  reaches  a specified  value  then  the  point 
at  the  fault  surface  enters  the  slip  routine.  Between  two 
consecutive  cycles,  n and  n+1,  the  accumulation  takes  the 
form 


n+1  n /IT  - 

e = e + - 


n + 1 n 
e = e 


) 


(3) 
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Rupture  occurs  if 


n+1 

c 


> W 


(4) 


where  W is  a specified  function  of  distance  from  the 
initial  point  of  rupture  (the  focus). 

Equation  (3)  is  similar  to  a plastic  work  criterion, 
where  the  plastic  work  (En+1)  is  given  by 


E 


n+1 


nn  Y2  /IT  - Y 

E Tv  ? 


(5) 


Equations  (3)  and  (5)  differ  only  by  the  factor  Y2/3p. 

We  have  been  successful  in  both  controlling  rupture 
velocity  and  reducing  the  stress  concentrations  at  the  end 
of  the  fault  by  allowing  W,  in  Eq.  (4),  to  be  a specified 
function  of  distance  from  the  point  of  rupture  initiation. 
The  functional  form  that  has  been  used  is 


= 6c  ) 

(x  + L/2  \2 

rj 

+ 

X 

f — i 
f-H 

l a—) 

L7 

ii 

4 ♦ s - 

tm 

+ L/2" 

“3 

0 < x 


x + 


(6a) 

(6b) 


where  L,  c and  d are  input  parameters.  The  rupture  is 
constrained  to  lie  between  - L/2  <_  x <_  L/2,  where  L is 
the  fault  length. 
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APPENDIX  D 

LINEAR  DIFFERENCING  SCHEME  USED  IN  LAGS 


The  following  discussion  concerns  the  differencing  scheme 
used  in  the  finite  difference  stress  wave  code  LAGS.  This  code 
performs  the  Linear  Analysis  of  geologic  Structures.  The  dif- 
ferencing is  similar  to  that  originally  presented  by  Wilkins^  and 
also  used  in  CRAM.  1 1 In  this  scheme  the  equilibrium  equations 
containing  the  stresses  and  their  spatial  partial  derivatives  are 
differenced.  The  stress -displacement  relationships  are  differenced 
separately  to  complete  the  finite-difference  formulation.  This 
procedure  is  superior  to  differencing  the  equilibrium  equations 
written  in  terns  of  displacements  for  two  important  reasons. 

First,  this  procedure  allows  the  treatment  of  almost  any  non- 
uniform  material  description.  Second,  boundary  condition  incorpora- 
tion is  facilitated  using  this  technique.  Both  axisymmetric  and 
plane  strain  configurations  can  be  considered. 

The  dynamic  equilibrium  equations  in  the  (x,y)  Cartesian 
coordinate  system  are  given  by 


pu  = 


3a, 

J 

3x” 


3a  a - a. 
X£  + x a 


3y 


pv  = 


3a  3a  a 

*£  ♦ __X  + JSZ 

3x  3y  x 


Cl) 


where  (u,v)  are  the  displacements  in  the  (x,y)  directions,  dots 
denote  time  derivatives,  p is  the  density,  6 is  the  circumferen- 
tial direction,  and  omn  is  the  stress  tensor.  The  * quantities 
are  non-zero  for  axisymmetric  configurations  only.  The  first  step 
in  the  determination  of  an  approximate  solution  of  Eq.  (1)  is  to 
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divide  the  continuum  into  a discrete  number  of  nodes.  As  depicted 
in  Fig.  1,  these  points  are  identified  in  the  two-dimensional  space 
by  their  i and  j values.  The  program  LAGS  assumes  a rectilinear 
mesh  geometry.  The  Dx^  and  Dy^  increments  describing  this 
geometry  (see  Fig.  1)  may  be  non-uniform  and  are  supplied  by  the 
user.  In  addition,  each  cell  in  the  mesh  must  be  supplied  a material 
identifier.  This  flexibility  allows  the  treatment  of  almost  any 
general  geologic  structure. 

If  Dt  is  the  time  increment  and  tn  a discrete  value  of 
time  (equal  to  nDt),  the  time  differenced  representation  of  Hq. 

(1)  is  based  on  the  following  analog: 

n+i  n-i 
a -a  + Dtiin 


an+1  = an  + Dta  ? , 

where  a is  u or  v.  The  spatial  difference  analogs  at  a parti- 
cular (i , j ) node  which  are  required  to  simulate  Eq.  (1)  are  of  the 
following  form:^*^ 

(?  ML5  of-  (v*.  - • V>\  * vO-  . 

t J » J 1 , J 

(3) 

fo  Iv)-  -=  7T — (a  Ax  + A Ax  - A Ax  - A Ax  \ 

\P  Pi,j  \ 1 1 3 3 4 4/ij 

where  the  subscript  n refers  to  the  nth  cell  of  Fig.  1 and  A 
is  some  component  of  stress.  The  terms  . , Axn,  and  Ayp  are 

defined  by: 


4>-  • 

i.J 

- £ 

k*l 

Ax 

* Ax 

= Dx . . 

i 

4 

l + l 

Ax 

* Ax 

= Dx. 
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in  a linear  isotropic  media  are  related  by 
'mn  '""kk^mn 


o__  - Ac,.,  <5  + 2^ie  , 

mn  * 


(6) 


"’here  A and  u are  the  Lame' constans , 6mn  is  the  Kronecker  delta 
function,  enn  is  the  strain  tensor,  and  double  subscripts  denote 
summation.  In  this  case 


e,,  = c + c + cn 
k k \ y Q 


The  linear  strain-displacement  relationships  are  given  by 


= iH 

"x  ~ ax 

_ av 
'y  ~ ay 

- 1 ( ~ + 3 v \ 

xy  I \3y  * ox) 

0,  plane  strain 


C7) 


e 


-,  axisymmetric 


Referring  to  Fig.  1,  the  strains  in  the  cell  bounded  by  the 
i*l>  i»  j - 1 , and  j-lines  are  given  by^>^ 


'x  “ TUx 


h (UM  ‘ Ui.M  • “i-1.3  • “i-l.J-l) 

r ‘ % fvi.J  * vi-i.i  ‘ vi,M  ‘ Vi.m) 

y ’ T _Ex~  (vi,j-l  * vi,j  • vi-i,j  ' vi-i , j - 1 ) 
* WJ  (ui,M  * “i-l.J-l  - “i.j  - Vij)] 


(8) 
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A>\  * = D>>; 


Ay  = Ay  = Dy . . 
3 •.  J 


For  axisymmetric  problems  the  °mn/px  terms  appearing  in  £4.  (1) 
are  represented  by: 


/A  \ _ 1 V Ak 

V p^/±  sr  ^77* 

,J  k=l  pkxk 


where  M represents  the  number  of  cells  surrounding  the  (i,j)  node 
and  Xj.  are  defined  as 


X = X = X . - 

1 * 1 + 1/2 


X = X = X.  . 

2 3 1-1/2 


Equations  (3)  through  (5)  ar  also  applicable  for  the  case 
when  (i , j ) is  not  an  interior  node.  To  illustrate  this  point,  the 
case  when  the  j-line  represents  a surface  in  space  where  a stress 
boundary  condition  is  applied  will  be  discussed.  Assuming  material 
exists  above  the  j-line,  cells  3 and  4 do  not  exist  for  this  case. 
Thus  the  summations  of  Eqs.  (4)  and  (3)  include  k = 1 and  k = 2 
only.  If  the  applied  pressure  is  P(x)  and  the  applied  shear 
stress  is  x(x)  Eq.  (3)  is  evaluated  using 


0>-3  = -p  (xi-l/2) 

°y»  * 'P  ^xi+l/2^ 

V.'  T Cxi-l/2^ 

°xy*'  1 <xiU/2>- 

Finite  dif ierencing  of  the  s tress -displacement  relationship 
completesthe  formulation.  From  Hooke's  law,  stresses  and  strains 
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and,  for  axisymnetric  geometry,  e is  calculated  using 


_ _ Av 

6 " v ex  ‘ e0 
o 


(9) 


where  the  bulk  strain  Av/v  is  evaluated  using  the  following 
expressions : 


Av  i T T + (x.  - l 
Av  _ ^ ! 2 \ i I 


DXi)  T3  * (\-  IDXi)' 


TT  ~ 


i i-1/2 


where 


T = Dx.Dy. 
i i 


T,  ’ 2ui,j-i  + 2ui-i,j  * ui(j  * “i-i.j.i 

T.  * D>i  (“i.j-l  - ui-l.j-l)  * Dxi  (vi-l,j  - 
T-  * % (Ui,j  ' ui-l,j)  * D*i  (vi,j  ' 

Equations  (2),  (3),  (4),  (5),  (8),  and  (9)  represent  the 
finite  difference  analog  of  Eqs.  (1),  (6),  and  (7). 
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AI PEND IX  E 


If  u,  v and  w are  the  x,  y and  z components  of  dis- 
placement, then  the  conservation  of  linear  momentum  may  be 
written 


3 CP  + cr  ) 9a 


u = 


p3x 


EL  ♦ 


3 a. 


xz 


pay 


p3z 


V 


3a  3 (P  + a ) 3a 

*y  + _ EL  * -ZL 


p J X 


poy 


pTz 


3a. 


w 


xz 


3a 


p3x  p3y 


XL  + 


3 (P  - a - a ) 

x y 

p'S  z 


(1) 

(2) 

(3) 


There  will  be  eight  elements  surrounding  a typical 
interior  nodal  point.  These  elements  are  labeled  A through 


H as  shown  in  Fig.  1.  The  exterior  points  of  element  A are 


labeled  1 through  8 as  shown 


Figure  1. 

Interior  points  of  8 ele- 
ments are  labeled  A,  B, 

C,  D,  E,  F.  G. 


in  Fig . 2 . 


Exterior  points  of  zone 
A are  labeled  1,  2,  3, 

4 , 5 , 6 , 7 , 8 . 
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The  following  calculations  are  performed  for  zone  A. 


l3x  I. 


(! 


A 

sy/A 


' 3 w \ 


p = js 
A F 


<°x> 


(V 

y A 


(°xy) 
xy  A 


(oxz> 


(oyz} 
yz  A 


u 

+ 

u 

+ u + u 

1 2 

1*  3 

5 6 6 7 

(4) 

Ax  + 

V 

+ 

V 

+ v + V 

1 4 

2 3 

6 7 6 8 

(5) 

Ay  + 

w 

+ 

w 

+ w + w 

1 s 

2 6 

3 7 4 8 

(6) 

Az  + 

/ 9u  \ 

/9v\  / 3w  \ 

(7) 

( 9x  ) 

4- 

\ 

(H  MJ 

y 

IT 

2 

(S) 

A ‘ Wa  ' 

(-)a] 

(8) 

y 

/ 9v 

/3u  \ 

(9) 

IT 

2 

L ‘ (3x'a  * i 

^)aJ 

u 

+ 

u + u + u 

0 

+ '0 

+ 0 +0 

^|oO 

II 

i_ 

u 

2 3 6 7 

5 8 + 

1 2 

•4  3 

5 6 8 7 

Ay+ 

Ax  + 

(10) 

u 

+ 

U + U + u 

w 

+ w 

+ w + w 

_ y 
8 

1 

5 

2 6 3 7 

JU.  4 

1 2 

4 3 

5 6 6 7 

Az  + 

Ax  + 

(ID 

V 

+ 

V + V + V 

W 

+ W 

+ w + w 

_ y 
" F 

1_ 

5 

2 6 3 7 

-Li  4 

1 4 

2 3 

6 7 5 8 

Az  + 

Ay  + 

(12) 


where  the  following  notation  has  been  adopted 


u =u-u  v = v - v , etc. 
12  1 2 12  1 2 


X - X 


y . y 


z - z 

3 7 


Ax+  = — - , Ay+  = -6— 2'  - 7 , Az  + = j 
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If  E is  a typical  stress  component  (P  + a , P + a 

^ " 7x  °y • aXy»  °x z * ayz^  ^ or  3»  then  the 

differenced  form  of  a typical  derivative  term  becomes 


1L  = ZAB  ^ 1 ' e7^1  ‘ eZ)  _ ZDCey(1  " eZ) 
p3x 


PAAX  + PBAX 


p^Ax  + pcAx 


Eef(1  - cy)ez 


r ,y  cz 
ahg  c e 


3 I 

p3y 


p3  z 


+ 

PpAx  + 

PpAx 

PHAx+ 

+ p^Ax 

. W1  - 

e ) (1  - 

e ) 

V X.  f m Z v 

"BCe  (!  ‘ e ) 

?aa/ 

* pDAy‘ 

pBAy+  + pCAy" 

. W1  • 

ex)  ez 

T 

FG 

eX  e 2 

pea/  - 

PHAy‘ 

pFAy+ 

* °GAy" 

_ rAE(1 

eX)(l  - 

£y) 

ZDHey  ^ 

°aaz* 

+ Pp  Az 

pdAz+  +phAz’ 

^ ZBFeX(^1 

+ 

- ey) 

+ 

T 

CG 

4. 

x y 

e e7 

Pp^6  P r-  A Z 


PCAZ  + pGAz 


(13) 


(14) 


(15) 


where 


<-x  = „y  _ Ay+  z Az  + 

b + - * e “ — + 7 * e " — t 

Ax  + Ax  Ay  + Ay  Az  + Az 

Equations  7 through  15  represent  finite  different  ap- 
proximations to  Hooke's  law  and  the  conservation  of  linear 
momentum.  In  order  to  write  the  equations  in  this  form, 
small  displacement  and  elastic  material  behavior  assumptions 
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required.  Nonlinear  material  behavior  will  be  simulated  in 
the  fault  zone  by  recasting  Eqs.  7 through  12  into  a stress 
rate-strain  rate  formulation.  This  will  entail  saving  six 
new  (stress)  variables  at  each  grid  location  lying  within 
the  fault  zone. 


